Machine learning project

Authors: Jose Pérez Cano & Álvaro Ribot Barrado

0. Libraries

install.packages("klaR")
install.packages("TunePareto")
install.packages("rgl")
install.packages("glmnet")
install.packages("ca")
# LDA/ QDA
library(MASS)

# RDA
library(klaR)

# Multinomial
library(nnet)

# Cross-Validation
library(TunePareto)

# Naive Bayes
library(e1071)

# k-NN
library(class)

# 3d
library(rgl)

# LASSO
library(Matrix)
library(glmnet)
Loading required package: foreach
Loaded glmnet 2.0-16
# Correspondence analysis
library(ca)

1. Read data

set.seed(2105)
setwd("../data")

The working directory was changed to /Users/joseperezcano/Desktop/CFIS segundo curso/AA1/Project/data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.
clev <- read.csv("cleveland.csv", header=F)

2. Preprocess data

We apply clustering and several plotting techniques to have an idea of the dataset. In case there are NAs we will use k-NN for imputation.

To extract new features we will apply PCA and FDA and keep this new features and components apart.

# It says which columns are all missings
# The index are returned in negative to eliminate them
na.columns <- function(dd){
  rmlist <- c()
  for (i in 1:ncol(dd)){
    if (min(dd[,i]) == max(dd[,i]) & min(dd[,i])==-9){
      rmlist <- c(rmlist, i)
    }
  }
  -rmlist
}

clev <- clev[,na.columns(clev)]

# Returns columns with more NA than a given threshold, also in negative
much.na.cols <- function(dd, threshold){
  rmlist <- c()
  for (i in 1:ncol(dd)){
    if (sum(dd[,i]==-9) > threshold){
      rmlist <- c(rmlist, i)
    }
  }   
  -rmlist  
}

clev <- clev[, much.na.cols(clev, 60)]

# Applies k-nearest neighbour imputation for a given variable
knn.imputation = function (dd, variable, varname, k)
{  
  aux = subset (dd, select = names(dd)[names(dd) != varname])
  aux1 = aux[!is.na(variable),]
  aux2 = aux[is.na(variable),]

  # Neither of aux1, aux2 can contain NAs
  knn.inc = knn (aux1,aux2, variable[!is.na(variable)], k)
  variable[is.na(variable)] = knn.inc
  variable
}

# This are the variables which values where substituted by dummy values.
dummy <- c("V1", "V2", "V36", "V69", "V70", "V71", "V72", "V73", "V28", "location")
clev <- clev[,!(names(clev) %in% dummy)]

# knn imputation for clev
na.names <- names(clev)[-much.na.cols(clev, 0)]
for (name in na.names){
  clev[, name][clev[, name] == -9] <- NA
  clev[, name] <- knn.imputation(clev, clev[,name], name, 7)
}

Now, we analyse the correlations among variables since it will have impact on later models.

corr.factors <- cor(clev)
which(abs(corr.factors)-diag(diag(corr.factors))>0.9, arr.ind=T)
    row col
V55  34  12
V57  36  14
V31  21  20
V29  20  21
V20  12  34
V22  14  36
rm.correlated <- c("V57", "V55")
clev <- clev[,!(names(clev) %in% rm.correlated)]

factores <- c("V58", "V4", "V9", "V16", "V18", "V19", "V20", "V21", "V22", "V23", "V24", "V25", "V26", "V27", "V38", "V39", "V41", "V51", "V56", "V11", "V59", "V60", "V61", "V63", "V65", "V67", "V68")
for (f in factores){
  clev[,f] <- as.factor(clev[,f])
}

# Dummy level gets replaced
clev$V25[clev$V25 == 2] <- 1
clev$V25 <- droplevels(clev$V25)
levels(clev$V25)
[1] "0" "1"

This is the dataset after the treatment for missings.

summary(clev)

       V3        V4      V9           V10        V11          V12       
 Min.   :29.00   0: 91   1: 22   Min.   : 94.0   0:108   Min.   :126.0  
 1st Qu.:48.00   1:191   2: 43   1st Qu.:120.0   1:174   1st Qu.:213.0  
 Median :55.00           3: 84   Median :130.0           Median :244.0  
 Mean   :54.41           4:133   Mean   :131.6           Mean   :249.1  
 3rd Qu.:61.00                   3rd Qu.:140.0           3rd Qu.:277.0  
 Max.   :77.00                   Max.   :200.0           Max.   :564.0  
                                                                        
      V14             V15        V16     V18     V19          V20    
 Min.   : 0.00   Min.   : 0.00   0:240   0:107   0:138   1      :38  
 1st Qu.: 0.00   1st Qu.: 0.00   1: 42   1:175   1:  2   12     :33  
 Median :10.00   Median :15.00                   2:142   2      :31  
 Mean   :16.64   Mean   :15.08                           8      :30  
 3rd Qu.:30.00   3rd Qu.:30.00                           6      :26  
 Max.   :99.00   Max.   :54.00                           11     :25  
                                                         (Other):99  
      V21      V22     V23     V24     V25     V26     V27          V29        
 4      : 13   81:67   0:271   0:186   0:211   0:252   0:248   Min.   : 1.800  
 15     : 13   82:96   1: 11   1: 95   1: 71   1: 30   1: 34   1st Qu.: 6.500  
 20     : 13   83:86           2:  1                           Median : 8.500  
 13     : 12   84:33                                           Mean   : 8.418  
 16     : 12                                                   3rd Qu.:10.075  
 21     : 12                                                   Max.   :15.000  
 (Other):207                                                                   
      V31              V32             V33              V34       
 Min.   : 3.000   Min.   : 71.0   Min.   : 40.00   Min.   : 84.0  
 1st Qu.: 7.000   1st Qu.:133.2   1st Qu.: 65.00   1st Qu.:154.0  
 Median : 9.000   Median :153.5   Median : 74.00   Median :168.0  
 Mean   : 9.754   Mean   :149.8   Mean   : 75.12   Mean   :168.1  
 3rd Qu.:12.000   3rd Qu.:165.8   3rd Qu.: 84.00   3rd Qu.:184.0  
 Max.   :18.000   Max.   :202.0   Max.   :119.00   Max.   :232.0  
                                                                  
      V35              V37         V38     V39          V40        V41    
 Min.   : 26.00   Min.   : 50.00   0:190   0:276   Min.   :0.000   1:135  
 1st Qu.: 70.00   1st Qu.: 80.00   1: 92   1:  6   1st Qu.:0.000   2:129  
 Median : 80.00   Median : 85.00                   Median :0.800   3: 18  
 Mean   : 78.74   Mean   : 84.95                   Mean   :1.027          
 3rd Qu.: 85.00   3rd Qu.: 90.00                   3rd Qu.:1.600          
 Max.   :120.00   Max.   :110.00                   Max.   :6.200          
                                                                          
      V43             V44         V51          V56      V58     V59     V60    
 Min.   : 24.0   Min.   :0.0000   1:  2   5      : 14   0:157   1:270   1:242  
 1st Qu.: 92.0   1st Qu.:0.0000   3:159   21     : 14   1: 50   2: 12   2: 40  
 Median :118.0   Median :0.0000   6: 14   14     : 12   2: 31                  
 Mean   :123.6   Mean   :0.6702   7:107   1      : 11   3: 32                  
 3rd Qu.:152.8   3rd Qu.:1.0000           17     : 11   4: 12                  
 Max.   :270.0   Max.   :3.0000           30     : 11                          
                                          (Other):209                          
 V61     V63     V65     V67     V68    
 1:224   1:238   1:236   1:233   1:246  
 2: 58   2: 44   2: 46   2: 49   2: 36

2.1 Visualizations

par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    hist(clev[,i], main = names(clev)[i], xlab="Values")
  }
}

par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    boxplot(clev[,i], xlab = names(clev)[i])
  }
}

par(mfrow = c(3,3))
for(i in 1:ncol(clev)){
  if (is.factor(clev[,i])) {
    hist(as.numeric(as.character(clev[,i])), main = names(clev)[i], xlab="Values")
  }
}

names_num <- c()
for(i in 1:ncol(clev)){
  if (!is.factor(clev[,i])) {
    names_num <- c(names_num, i)
  }
}
clev_numeric <- clev[,names_num]

clev_cor <- cor(clev_numeric)
which(clev_cor > 0.5 & clev_cor < 1, arr.ind = TRUE)
    row col
V34  10   2
V37  12   2
V15   5   4
V14   4   5
V31   7   6
V29   6   7
V32   8   7
V31   7   8
V10   2  10
V37  12  11
V10   2  12
V35  11  12
which(-clev_cor > 0.5 & -clev_cor < 1, arr.ind = TRUE)
     row col

2.2 Modification of values

par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
    qqnorm(clev_numeric[,i], main = c("Q-Q Plot: ", names(clev_numeric)[i]))
    qqline(clev_numeric[,i], col=2)
}

Boxcox

par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
  boxcox(lm(clev_numeric[,i]-min(clev_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(clev_numeric))[i])
}

par(mfrow = c(1,3))
#we treat them as special cases
aux <- clev_numeric[,"V14"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
aux <- clev_numeric[,"V15"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
aux <- clev_numeric[,"V40"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))

clev_sqrt <- c("V10", "V12", "V31", "V43")
clev_sqrt_especial <- c("V14", "V40")
clev_box <- clev_numeric
#box-cox transformation
for (i in 1:ncol(clev_box)){
  if (names(clev_box)[i] %in% clev_sqrt) {
    clev_box[,i] <- 2*sqrt(clev_box[,i]-min(clev_box[,i])+1e-6)
  } else if (names(clev_box)[i] %in% clev_sqrt_especial){
    clev_box[,i] <- 2*sqrt(clev_box[,i])
  }
}
aux <- names(clev_box)
clev_box <- data.frame(clev_box) 
colnames <- aux
par(mfrow = c(2,3))
for(i in 1:ncol(clev_box)){
    qqnorm(clev_box[,i], main = c("Q-Q Plot: ", names(clev_box)[i]))
    qqline(clev_box[,i], col=2)
}

For further models, normalisation will be useful.

And this is the final dataset after processing it.

clev[, colnames(clev_box)] <- clev_box[,colnames(clev_box)]
summary(clev)
       V3           V4      V9           V10           V11    
 Min.   :-2.80693   0: 91   1: 22   Min.   :-3.90175   0:108  
 1st Qu.:-0.70820   1:191   2: 43   1st Qu.:-0.55123   1:174  
 Median : 0.06502           3: 84   Median : 0.04091          
 Mean   : 0.00000           4:133   Mean   : 0.00000          
 3rd Qu.: 0.72778                   3rd Qu.: 0.55507          
 Max.   : 2.49513                   Max.   : 2.86408          
                                                              
      V12                 V14               V15           V16     V18     V19    
 Min.   :-4.643087   Min.   :-1.0513   Min.   :-0.98383   0:240   0:107   0:138  
 1st Qu.:-0.650497   1st Qu.:-1.0513   1st Qu.:-0.98383   1: 42   1:175   1:  2  
 Median : 0.006803   Median : 0.0726   Median :-0.00532                   2:142  
 Mean   : 0.000000   Mean   : 0.0000   Mean   : 0.00000                          
 3rd Qu.: 0.617022   3rd Qu.: 0.8954   3rd Qu.: 0.97319                          
 Max.   : 4.315880   Max.   : 2.4850   Max.   : 2.53880                          
                                                                                 
      V20          V21      V22     V23     V24     V25     V26     V27    
 1      :38   4      : 13   81:67   0:271   0:186   0:211   0:252   0:248  
 12     :33   15     : 13   82:96   1: 11   1: 95   1: 71   1: 30   1: 34  
 2      :31   20     : 13   83:86           2:  1                          
 8      :30   13     : 12   84:33                                          
 6      :26   16     : 12                                                  
 11     :25   21     : 12                                                  
 (Other):99   (Other):207                                                  
      V29                V31               V32               V33          
 Min.   :-2.55482   Min.   :-4.1234   Min.   :-3.4360   Min.   :-2.54549  
 1st Qu.:-0.74055   1st Qu.:-0.8588   1st Qu.:-0.7205   1st Qu.:-0.73334  
 Median : 0.03148   Median :-0.1247   Median : 0.1629   Median :-0.08097  
 Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
 3rd Qu.: 0.63946   3rd Qu.: 0.7743   3rd Qu.: 0.6973   3rd Qu.: 0.64389  
 Max.   : 2.54059   Max.   : 2.2000   Max.   : 2.2786   Max.   : 3.18089  
                                                                          
      V34                 V35                V37            V38     V39    
 Min.   :-3.565358   Min.   :-3.95648   Min.   :-3.687203   0:190   0:276  
 1st Qu.:-0.596232   1st Qu.:-0.65596   1st Qu.:-0.521933   1: 92   1:  6  
 Median :-0.002407   Median : 0.09416   Median : 0.005612                  
 Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.000000                  
 3rd Qu.: 0.676251   3rd Qu.: 0.46922   3rd Qu.: 0.533157                  
 Max.   : 2.712223   Max.   : 3.09464   Max.   : 2.643337                  
                                                                           
      V40          V41          V43                V44          V51    
 Min.   :-1.1986   1:135   Min.   :-3.73286   Min.   :-0.7158   1:  2  
 1st Qu.:-1.1986   2:129   1st Qu.:-0.53995   1st Qu.:-0.7158   3:159  
 Median : 0.1781   3: 18   Median : 0.02122   Median :-0.7158   6: 14  
 Mean   : 0.0000           Mean   : 0.00000   Mean   : 0.0000   7:107  
 3rd Qu.: 0.7484           3rd Qu.: 0.66073   3rd Qu.: 0.3522          
 Max.   : 2.6341           Max.   : 2.34044   Max.   : 2.4883          
                                                                       
      V56      V58     V59     V60     V61     V63     V65     V67     V68    
 5      : 14   0:157   1:270   1:242   1:224   1:238   1:236   1:233   1:246  
 21     : 14   1: 50   2: 12   2: 40   2: 58   2: 44   2: 46   2: 49   2: 36  
 14     : 12   2: 31                                                          
 1      : 11   3: 32                                                          
 17     : 11   4: 12                                                          
 30     : 11                                                                  
 (Other):209

2.3. Feature extraction

Separe train and test data, seed for reproducibility.

set.seed(2000)
n <- nrow(clev)
train.lenght <- round(2*n/3)

clev <- clev[sample(n),]
train <- clev[1:train.lenght,]
test <- clev[(train.lenght+1):n,]
names_num <- c()
for(i in 1:ncol(train)){
  if (!is.factor(train[,i])) {
    names_num <- c(names_num, i)
  }
}
train_num <- train[,names_num]

Extract PCA features.

pca <- princomp(train_num)
screeplot(pca)

summary(pca)
Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4
Standard deviation     1.8444732 1.5577035 1.3875402 1.21213461
Proportion of Variance 0.2297381 0.1638543 0.1300108 0.09921788
Cumulative Proportion  0.2297381 0.3935924 0.5236032 0.62282104
                           Comp.5     Comp.6     Comp.7    Comp.8
Standard deviation     1.01035405 0.97964763 0.85803949 0.7938423
Proportion of Variance 0.06893431 0.06480791 0.04971676 0.0425556
Cumulative Proportion  0.69175535 0.75656326 0.80628002 0.8488356
                           Comp.9    Comp.10    Comp.11
Standard deviation     0.73133601 0.69242679 0.62465621
Proportion of Variance 0.03611787 0.03237695 0.02634938
Cumulative Proportion  0.88495350 0.91733045 0.94367983
                          Comp.12    Comp.13    Comp.14
Standard deviation     0.57654420 0.50756673 0.43609672
Proportion of Variance 0.02244675 0.01739701 0.01284263
Cumulative Proportion  0.96612658 0.98352359 0.99636622
                           Comp.15
Standard deviation     0.231971913
Proportion of Variance 0.003633784
Cumulative Proportion  1.000000000

Fp <- pca$scores
Gs <- pca$loadings

Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) * 2

col.class <- as.numeric(train$V58)
col.class[col.class==1] <- "red"
col.class[col.class==2] <- "green"
col.class[col.class==3] <- "blue"
col.class[col.class==4] <- "yellow"
col.class[col.class==5] <- "purple"

plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])
text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))

V1, V59, V57 creates many problems

problematic <- c("V57", "V59")
train <- train[, !(names(train) %in% problematic)]
fda <- lda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))

train$LD1 <- loadings[,1]
train$LD2 <- loadings[,2]
train$LD3 <- loadings[,3]
train$LD4 <- loadings[,4]

fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]
test$LD2 <- fda_test$x[,2]
test$LD3 <- fda_test$x[,3]
test$LD4 <- fda_test$x[,4]

Análisis por correspondencias

ac <- mjca(clev[,names(clev) %in% factores], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")


ac_ind <- mjca(train[,names(train) %in% factores], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))


mca.features <- ac_ind$rowcoord

3. Resampling protocol

Principal utility function for doing cross-validation.

Create the CV function for any model with associated predict and update functions.

cross.validation <- function(data, target, model, times, nfolds, need.class){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      print(paste0("Fold: ", j))
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      model <- update(model, data=tr)
      pred <- predict(model, newdata=va)
      if (need.class){
        pred <- pred$class
      }
      
      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
    print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
  }
  mean(err.total)
}

Cross-validation for k-Nearest Neighbour

cross.validation.knn <- function(data, target, times, nfolds, K){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      pred <- knn(tr, va, target[-val], k = K)

      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
  }
  mean(err.total)
}

Cross-validation Naive-Bayes

cross.validation.naive <- function(data, target, model, times, nfolds){
  set.seed(0202)
  CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)

  err.total <- c()
  for (i in 1:times){
    err.onetime <- c()
    for (j in 1:nfolds){
      print(paste0("Fold: ", j))
      val <- unlist(CV.folds[[i]][[j]])
      
      tr <- data[-val,]
      va <- data[val,]

      model <- naiveBayes(V58~.-LD1-LD2-LD3-LD4, data=tr)
      pred <- predict(model, newdata=va)
      
      err.table <- table(True=target[val], Predicted=pred)
      err <- 1-sum(diag(err.table))/sum(err.table)
      err.onetime <- c(err.onetime, err)
    }
    err.total <- c(err.total, mean(err.onetime))
    print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
  }
  mean(err.total)
}

4. Models

The models we are going to use are: - LDA - QDA - RDA - k-NN - Naïve Bayes - GLM

rda.model <- rda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
naive.model <- naiveBayes(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
cross.validation(train, train$V58, rda.model, 10, 10, T)
rda.model.fda <- rda(V58~.,data=train)
cross.validation(train, train$V58, rda.model.fda, 10, 10, T)
cross.validation.naive(train, train$V58, naive.model, 10, 10)
err <- c()
for (k in 1:20){
  err <- c(err, cross.validation.knn(train, train$V58, 10,10, k))
}
plot(err, type = "l")

err
 [1] 0.1519883 0.1965205 0.1950585 0.2004678 0.2073977 0.2152924
 [7] 0.2184795 0.2276023 0.2313450 0.2323392 0.2419883 0.2456433
[13] 0.2642398 0.2690351 0.2743275 0.2791228 0.2871053 0.2955556
[19] 0.3003801 0.3050877
cross.validation.knn(train, train$V58, 10, 10, 1)

multinomial.model <- multinom(V58~., data=train)

cross.validation(train, train$V58, multinomial.model, 10, 10, F)

multinomial.model.step <- step(multinomial.model)

cross.validation(train, train$V58, multinomial.model.step, 10, 10, F)

multinomial.model.noFDA <- multinom(V58~.-LD1-LD2-LD3-LD4, data=train)

cross.validation(train, train$V58, multinomial.model.noFDA, 10, 10, F)

multinomial.model.noFDA.step <- step(multinomial.model.noFDA)

cross.validation(train, train$V58, multinomial.model.noFDA.step, 10, 10, F)

Test error

rda.model <- update(rda.model.fda, data=train)
pred.test <- predict(rda.model.fda, test)
pred.test <- pred.test$class
(err.table <- table(True=test$V58, Pred=pred.test))
    Pred
True  0  1  2  3  4
   0 40  5  0  0  0
   1  6 13  1  0  2
   2  0  2  6  4  1
   3  0  0  1  9  0
   4  0  0  0  1  3
(err.test <- 1-sum(diag(err.table))/sum(err.table))
[1] 0.2446809

Parte 2

We will now use model on a dataset not yet investigated: Hungary.

hung <- read.table("../data/processed.hungarian.data", header=F, sep=',', na.strings="?")
summary(hung)
       V1              V2               V3              V4       
 Min.   :28.00   Min.   :0.0000   Min.   :1.000   Min.   : 92.0  
 1st Qu.:42.00   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:120.0  
 Median :49.00   Median :1.0000   Median :3.000   Median :130.0  
 Mean   :47.83   Mean   :0.7245   Mean   :2.983   Mean   :132.6  
 3rd Qu.:54.00   3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:140.0  
 Max.   :66.00   Max.   :1.0000   Max.   :4.000   Max.   :200.0  
                                                  NA's   :1      
       V5              V6                V7        
 Min.   : 85.0   Min.   :0.00000   Min.   :0.0000  
 1st Qu.:209.0   1st Qu.:0.00000   1st Qu.:0.0000  
 Median :243.0   Median :0.00000   Median :0.0000  
 Mean   :250.8   Mean   :0.06993   Mean   :0.2184  
 3rd Qu.:282.5   3rd Qu.:0.00000   3rd Qu.:0.0000  
 Max.   :603.0   Max.   :1.00000   Max.   :2.0000  
 NA's   :23      NA's   :8         NA's   :1       
       V8              V9              V10        
 Min.   : 82.0   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:122.0   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :140.0   Median :0.0000   Median :0.0000  
 Mean   :139.1   Mean   :0.3038   Mean   :0.5861  
 3rd Qu.:155.0   3rd Qu.:1.0000   3rd Qu.:1.0000  
 Max.   :190.0   Max.   :1.0000   Max.   :5.0000  
 NA's   :1       NA's   :1                        
      V11             V12           V13             V14        
 Min.   :1.000   Min.   :0     Min.   :3.000   Min.   :0.0000  
 1st Qu.:2.000   1st Qu.:0     1st Qu.:5.250   1st Qu.:0.0000  
 Median :2.000   Median :0     Median :6.000   Median :0.0000  
 Mean   :1.894   Mean   :0     Mean   :5.643   Mean   :0.3605  
 3rd Qu.:2.000   3rd Qu.:0     3rd Qu.:7.000   3rd Qu.:1.0000  
 Max.   :3.000   Max.   :0     Max.   :7.000   Max.   :1.0000  
 NA's   :190     NA's   :291   NA's   :266

Missing values treatment

rm.var <- c("V11", "V12", "V13")
hung <- hung[,!(names(hung) %in% rm.var)]
hung[is.na(hung)] <- -9
var.imputable <- c("V4", "V5", "V6", "V7", "V8", "V9")
for (nom in var.imputable){
  hung[, nom][hung[, nom] == -9] <- NA
  variable <- hung[, nom]
  hung[, nom] <- knn.imputation(hung, variable, nom, 7)
}
var.factor <- c("V2", "V3", "V6", "V7", "V9", "V14")
for (nom in var.factor){
  hung[, nom] <- as.factor(as.character(hung[,nom]))
}
summary(hung)

       V1        V2      V3            V4              V5       
 Min.   :28.00   0: 81   1: 11   Min.   : 27.0   Min.   :  4.0  
 1st Qu.:42.00   1:213   2:106   1st Qu.:120.0   1st Qu.:198.0  
 Median :49.00           3: 54   Median :130.0   Median :237.0  
 Mean   :47.83           4:123   Mean   :132.2   Mean   :236.6  
 3rd Qu.:54.00                   3rd Qu.:140.0   3rd Qu.:277.0  
 Max.   :66.00                   Max.   :200.0   Max.   :603.0  
 V6      V7            V8        V9           V10         V14    
 0:266   0:235   Min.   : 54.0   0:204   Min.   :0.0000   0:188  
 1: 28   1: 53   1st Qu.:122.0   1: 89   1st Qu.:0.0000   1:106  
         2:  6   Median :140.0   2:  1   Median :0.0000          
                 Mean   :138.8           Mean   :0.5861          
                 3rd Qu.:155.0           3rd Qu.:1.0000          
                 Max.   :190.0           Max.   :5.0000

Gaussianity treatment

names_num <- c()
for(i in 1:ncol(hung)){
  if (!is.factor(hung[,i])) {
    names_num <- c(names_num, i)
  }
}
hung_numeric <- hung[,names_num]

hung_cor <- cor(hung_numeric)
which(hung_cor > 0.5 & hung_cor < 1, arr.ind = TRUE)
     row col
which(-hung_cor > 0.5 & -hung_cor < 1, arr.ind = TRUE)
     row col
par(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
    qqnorm(hung_numeric[,i], main = c("Q-Q Plot: ", names(hung_numeric)[i]))
    qqline(hung_numeric[,i], col=2)
}

Boxcox

par(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
  boxcox(lm(hung_numeric[,i]-min(hung_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(hung_numeric))[i])
}

We see var 10 requires a logaritmic transformation.

hung[,"V10"] <- log(1+hung[,"V10"])
qqnorm(hung[,"V10"])
qqline(hung[,"V10"], col=2)

Feature extraction

Separe train and test data, seed for reproducibility.

set.seed(2000)
n <- nrow(hung)
train.lenght <- round(2*n/3)

hung <- hung[sample(n),]
train <- hung[1:train.lenght,]
test <- hung[(train.lenght+1):n,]
names_num <- c()
for(i in 1:ncol(train)){
  if (!is.factor(train[,i])) {
    names_num <- c(names_num, i)
  }
}
train_num <- train[,names_num]

Extract PCA features.

pca <- princomp(train_num)
screeplot(pca)

summary(pca)
Importance of components:
                           Comp.1      Comp.2      Comp.3
Standard deviation     86.0161788 24.47912423 17.80405417
Proportion of Variance  0.8849996  0.07167612  0.03791583
Cumulative Proportion   0.8849996  0.95667568  0.99459151
                            Comp.4       Comp.5
Standard deviation     6.709559731 0.4448637716
Proportion of Variance 0.005384815 0.0000236721
Cumulative Proportion  0.999976328 1.0000000000

Fp <- pca$scores
Gs <- pca$loadings

Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) *0.3

col.class <- as.numeric(train$V14)
col.class[col.class==0] <- "red"
col.class[col.class==1] <- "blue"

plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])
text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))

fda <- lda(V14~., data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))

train$LD1 <- loadings[,1]

fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]

Análisis por correspondencias

ac <- mjca(hung[,names(hung) %in% var.factor], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")


ac_ind <- mjca(train[,names(train) %in% var.factor], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))


mca.features <- ac_ind$rowcoord

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCiMgTWFjaGluZSBsZWFybmluZyBwcm9qZWN0CiMjIEF1dGhvcnM6IEpvc2UgUMOpcmV6IENhbm8gJiDDgWx2YXJvIFJpYm90IEJhcnJhZG8KCiMjIyAwLiBMaWJyYXJpZXMKCmBgYHtyfQppbnN0YWxsLnBhY2thZ2VzKCJrbGFSIikKaW5zdGFsbC5wYWNrYWdlcygiVHVuZVBhcmV0byIpCmluc3RhbGwucGFja2FnZXMoInJnbCIpCmluc3RhbGwucGFja2FnZXMoImdsbW5ldCIpCmluc3RhbGwucGFja2FnZXMoImNhIikKYGBgCgpgYGB7cn0KIyBMREEvIFFEQQpsaWJyYXJ5KE1BU1MpCgojIFJEQQpsaWJyYXJ5KGtsYVIpCgojIE11bHRpbm9taWFsCmxpYnJhcnkobm5ldCkKCiMgQ3Jvc3MtVmFsaWRhdGlvbgpsaWJyYXJ5KFR1bmVQYXJldG8pCgojIE5haXZlIEJheWVzCmxpYnJhcnkoZTEwNzEpCgojIGstTk4KbGlicmFyeShjbGFzcykKCiMgM2QKbGlicmFyeShyZ2wpCgojIExBU1NPCmxpYnJhcnkoTWF0cml4KQpsaWJyYXJ5KGdsbW5ldCkKCiMgQ29ycmVzcG9uZGVuY2UgYW5hbHlzaXMKbGlicmFyeShjYSkKYGBgCgoKIyMjIDEuIFJlYWQgZGF0YQoKYGBge3J9CnNldC5zZWVkKDIxMDUpCnNldHdkKCIuLi9kYXRhIikKY2xldiA8LSByZWFkLmNzdigiY2xldmVsYW5kLmNzdiIsIGhlYWRlcj1GKQpodW5nIDwtIHJlYWQuY3N2KCJodW5nYXJpYW4uY3N2IiwgaGVhZGVyPUYpCnZhIDwtIHJlYWQuY3N2KCJsb25nLWJlYWNoLXZhLmNzdiIsIGhlYWRlcj1GKQpzd2l0eiA8LSByZWFkLmNzdigic3dpdHplcmxhbmQuY3N2IiwgaGVhZGVyPUYpCgpjbGV2JGxvY2F0aW9uIDwtICJjbGV2ZWxhbmQiCmh1bmckbG9jYXRpb24gPC0gImh1bmdhcmlhbiIKdmEkbG9jYXRpb24gPC0gImxvbmctYmVhY2gtdmEiCnN3aXR6JGxvY2F0aW9uIDwtICJzd2l0emVybGFuZCIKCmhlYXJ0MSA8LSByYmluZChjbGV2LCBodW5nKQpoZWFydDIgPC0gcmJpbmQodmEsIHN3aXR6KQpoZWFydCA8LSByYmluZChoZWFydDEsIGhlYXJ0MikKaGVhZChoZWFydCkKYGBgCgoKIyMjIDIuIFByZXByb2Nlc3MgZGF0YQpXZSBhcHBseSBjbHVzdGVyaW5nIGFuZCBzZXZlcmFsIHBsb3R0aW5nIHRlY2huaXF1ZXMgdG8gaGF2ZSBhbiBpZGVhIG9mIHRoZSBkYXRhc2V0LiBJbiBjYXNlIHRoZXJlIGFyZSBOQXMgd2Ugd2lsbCB1c2Ugay1OTiBmb3IgaW1wdXRhdGlvbi4gCgpUbyBleHRyYWN0IG5ldyBmZWF0dXJlcyB3ZSB3aWxsIGFwcGx5IFBDQSBhbmQgRkRBIGFuZCBrZWVwIHRoaXMgbmV3IGZlYXR1cmVzIGFuZCBjb21wb25lbnRzIGFwYXJ0LgoKYGBge3J9CiMgSXQgc2F5cyB3aGljaCBjb2x1bW5zIGFyZSBhbGwgbWlzc2luZ3MKIyBUaGUgaW5kZXggYXJlIHJldHVybmVkIGluIG5lZ2F0aXZlIHRvIGVsaW1pbmF0ZSB0aGVtCm5hLmNvbHVtbnMgPC0gZnVuY3Rpb24oZGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAobWluKGRkWyxpXSkgPT0gbWF4KGRkWyxpXSkgJiBtaW4oZGRbLGldKT09LTkpewogICAgICBybWxpc3QgPC0gYyhybWxpc3QsIGkpCiAgICB9CiAgfQogIC1ybWxpc3QKfQoKY2xldiA8LSBjbGV2WyxuYS5jb2x1bW5zKGNsZXYpXQoKIyBSZXR1cm5zIGNvbHVtbnMgd2l0aCBtb3JlIE5BIHRoYW4gYSBnaXZlbiB0aHJlc2hvbGQsIGFsc28gaW4gbmVnYXRpdmUKbXVjaC5uYS5jb2xzIDwtIGZ1bmN0aW9uKGRkLCB0aHJlc2hvbGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAoc3VtKGRkWyxpXT09LTkpID4gdGhyZXNob2xkKXsKICAgICAgcm1saXN0IDwtIGMocm1saXN0LCBpKQogICAgfQogIH0gICAKICAtcm1saXN0ICAKfQoKY2xldiA8LSBjbGV2WywgbXVjaC5uYS5jb2xzKGNsZXYsIDYwKV0KCiMgQXBwbGllcyBrLW5lYXJlc3QgbmVpZ2hib3VyIGltcHV0YXRpb24gZm9yIGEgZ2l2ZW4gdmFyaWFibGUKa25uLmltcHV0YXRpb24gPSBmdW5jdGlvbiAoZGQsIHZhcmlhYmxlLCB2YXJuYW1lLCBrKQp7ICAKICBhdXggPSBzdWJzZXQgKGRkLCBzZWxlY3QgPSBuYW1lcyhkZClbbmFtZXMoZGQpICE9IHZhcm5hbWVdKQogIGF1eDEgPSBhdXhbIWlzLm5hKHZhcmlhYmxlKSxdCiAgYXV4MiA9IGF1eFtpcy5uYSh2YXJpYWJsZSksXQoKICAjIE5laXRoZXIgb2YgYXV4MSwgYXV4MiBjYW4gY29udGFpbiBOQXMKICBrbm4uaW5jID0ga25uIChhdXgxLGF1eDIsIHZhcmlhYmxlWyFpcy5uYSh2YXJpYWJsZSldLCBrKQogIHZhcmlhYmxlW2lzLm5hKHZhcmlhYmxlKV0gPSBrbm4uaW5jCiAgdmFyaWFibGUKfQoKIyBUaGlzIGFyZSB0aGUgdmFyaWFibGVzIHdoaWNoIHZhbHVlcyB3aGVyZSBzdWJzdGl0dXRlZCBieSBkdW1teSB2YWx1ZXMuCmR1bW15IDwtIGMoIlYxIiwgIlYyIiwgIlYzNiIsICJWNjkiLCAiVjcwIiwgIlY3MSIsICJWNzIiLCAiVjczIiwgIlYyOCIsICJsb2NhdGlvbiIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIGR1bW15KV0KCiMga25uIGltcHV0YXRpb24gZm9yIGNsZXYKbmEubmFtZXMgPC0gbmFtZXMoY2xldilbLW11Y2gubmEuY29scyhjbGV2LCAwKV0KZm9yIChuYW1lIGluIG5hLm5hbWVzKXsKICBjbGV2WywgbmFtZV1bY2xldlssIG5hbWVdID09IC05XSA8LSBOQQogIGNsZXZbLCBuYW1lXSA8LSBrbm4uaW1wdXRhdGlvbihjbGV2LCBjbGV2WyxuYW1lXSwgbmFtZSwgNykKfQpgYGAKCk5vdywgd2UgYW5hbHlzZSB0aGUgY29ycmVsYXRpb25zIGFtb25nIHZhcmlhYmxlcyBzaW5jZSBpdCB3aWxsIGhhdmUgaW1wYWN0IG9uIGxhdGVyIG1vZGVscy4KCmBgYHtyfQpjb3JyLmZhY3RvcnMgPC0gY29yKGNsZXYpCndoaWNoKGFicyhjb3JyLmZhY3RvcnMpLWRpYWcoZGlhZyhjb3JyLmZhY3RvcnMpKT4wLjksIGFyci5pbmQ9VCkKCnJtLmNvcnJlbGF0ZWQgPC0gYygiVjU3IiwgIlY1NSIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIHJtLmNvcnJlbGF0ZWQpXQoKZmFjdG9yZXMgPC0gYygiVjU4IiwgIlY0IiwgIlY5IiwgIlYxNiIsICJWMTgiLCAiVjE5IiwgIlYyMCIsICJWMjEiLCAiVjIyIiwgIlYyMyIsICJWMjQiLCAiVjI1IiwgIlYyNiIsICJWMjciLCAiVjM4IiwgIlYzOSIsICJWNDEiLCAiVjUxIiwgIlY1NiIsICJWMTEiLCAiVjU5IiwgIlY2MCIsICJWNjEiLCAiVjYzIiwgIlY2NSIsICJWNjciLCAiVjY4IikKZm9yIChmIGluIGZhY3RvcmVzKXsKICBjbGV2WyxmXSA8LSBhcy5mYWN0b3IoY2xldlssZl0pCn0KCiMgRHVtbXkgbGV2ZWwgZ2V0cyByZXBsYWNlZApjbGV2JFYyNVtjbGV2JFYyNSA9PSAyXSA8LSAxCmNsZXYkVjI1IDwtIGRyb3BsZXZlbHMoY2xldiRWMjUpCmxldmVscyhjbGV2JFYyNSkKYGBgCgpUaGlzIGlzIHRoZSBkYXRhc2V0IGFmdGVyIHRoZSB0cmVhdG1lbnQgZm9yIG1pc3NpbmdzLgoKYGBge3J9CnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjEgVmlzdWFsaXphdGlvbnMKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBoaXN0KGNsZXZbLGldLCBtYWluID0gbmFtZXMoY2xldilbaV0sIHhsYWI9IlZhbHVlcyIpCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBib3hwbG90KGNsZXZbLGldLCB4bGFiID0gbmFtZXMoY2xldilbaV0pCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDMsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKGlzLmZhY3RvcihjbGV2WyxpXSkpIHsKICAgIGhpc3QoYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoY2xldlssaV0pKSwgbWFpbiA9IG5hbWVzKGNsZXYpW2ldLCB4bGFiPSJWYWx1ZXMiKQogIH0KfQpgYGAKCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2woY2xldikpewogIGlmICghaXMuZmFjdG9yKGNsZXZbLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQpjbGV2X251bWVyaWMgPC0gY2xldlssbmFtZXNfbnVtXQoKY2xldl9jb3IgPC0gY29yKGNsZXZfbnVtZXJpYykKd2hpY2goY2xldl9jb3IgPiAwLjUgJiBjbGV2X2NvciA8IDEsIGFyci5pbmQgPSBUUlVFKQp3aGljaCgtY2xldl9jb3IgPiAwLjUgJiAtY2xldl9jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKYGBgCgoKIyMjIyAyLjIgTW9kaWZpY2F0aW9uIG9mIHZhbHVlcwoKYGBge3J9CnBhcihtZnJvdyA9IGMoMiwzKSkKZm9yKGkgaW4gMTpsZW5ndGgoY2xldl9udW1lcmljKSl7CiAgICBxcW5vcm0oY2xldl9udW1lcmljWyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X251bWVyaWMpW2ldKSkKICAgIHFxbGluZShjbGV2X251bWVyaWNbLGldLCBjb2w9MikKfQpgYGAKCgpCb3hjb3gKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGNsZXZfbnVtZXJpYykpewogIGJveGNveChsbShjbGV2X251bWVyaWNbLGldLW1pbihjbGV2X251bWVyaWNbLGldKSsxZS02fjEpLGxhbWJkYSA9IHNlcSgtMSwgMS41LCBieT0wLjEpLCB4bGFiID0gYyhuYW1lcyhjbGV2X251bWVyaWMpKVtpXSkKfQpgYGAKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDEsMykpCiN3ZSB0cmVhdCB0aGVtIGFzIHNwZWNpYWwgY2FzZXMKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNSJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlY0MCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYGBgCgoKYGBge3J9CmNsZXZfc3FydCA8LSBjKCJWMTAiLCAiVjEyIiwgIlYzMSIsICJWNDMiKQpjbGV2X3NxcnRfZXNwZWNpYWwgPC0gYygiVjE0IiwgIlY0MCIpCmNsZXZfYm94IDwtIGNsZXZfbnVtZXJpYwojYm94LWNveCB0cmFuc2Zvcm1hdGlvbgpmb3IgKGkgaW4gMTpuY29sKGNsZXZfYm94KSl7CiAgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydCkgewogICAgY2xldl9ib3hbLGldIDwtIDIqc3FydChjbGV2X2JveFssaV0tbWluKGNsZXZfYm94WyxpXSkrMWUtNikKICB9IGVsc2UgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydF9lc3BlY2lhbCl7CiAgICBjbGV2X2JveFssaV0gPC0gMipzcXJ0KGNsZXZfYm94WyxpXSkKICB9Cn0KYXV4IDwtIG5hbWVzKGNsZXZfYm94KQpjbGV2X2JveCA8LSBkYXRhLmZyYW1lKGNsZXZfYm94KSAKY29sbmFtZXMgPC0gYXV4CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2X2JveCkpewogICAgcXFub3JtKGNsZXZfYm94WyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X2JveClbaV0pKQogICAgcXFsaW5lKGNsZXZfYm94WyxpXSwgY29sPTIpCn0KYGBgCgpGb3IgZnVydGhlciBtb2RlbHMsIG5vcm1hbGlzYXRpb24gd2lsbCBiZSB1c2VmdWwuCgpgYGB7cn0KY2xldl9ib3ggPC0gYXBwbHkoY2xldl9ib3gsIDIsIHNjYWxlKQpgYGAKCgpBbmQgdGhpcyBpcyB0aGUgZmluYWwgZGF0YXNldCBhZnRlciBwcm9jZXNzaW5nIGl0LgoKYGBge3J9CmNsZXZbLCBjb2xuYW1lcyhjbGV2X2JveCldIDwtIGNsZXZfYm94Wyxjb2xuYW1lcyhjbGV2X2JveCldCnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjMuIEZlYXR1cmUgZXh0cmFjdGlvbgoKU2VwYXJlIHRyYWluIGFuZCB0ZXN0IGRhdGEsIHNlZWQgZm9yIHJlcHJvZHVjaWJpbGl0eS4KCmBgYHtyfQpzZXQuc2VlZCgyMDAwKQpuIDwtIG5yb3coY2xldikKdHJhaW4ubGVuZ2h0IDwtIHJvdW5kKDIqbi8zKQoKY2xldiA8LSBjbGV2W3NhbXBsZShuKSxdCnRyYWluIDwtIGNsZXZbMTp0cmFpbi5sZW5naHQsXQp0ZXN0IDwtIGNsZXZbKHRyYWluLmxlbmdodCsxKTpuLF0KYGBgCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2wodHJhaW4pKXsKICBpZiAoIWlzLmZhY3Rvcih0cmFpblssaV0pKSB7CiAgICBuYW1lc19udW0gPC0gYyhuYW1lc19udW0sIGkpCiAgfQp9CnRyYWluX251bSA8LSB0cmFpblssbmFtZXNfbnVtXQpgYGAKCkV4dHJhY3QgUENBIGZlYXR1cmVzLgoKYGBge3J9CnBjYSA8LSBwcmluY29tcCh0cmFpbl9udW0pCnNjcmVlcGxvdChwY2EpCnN1bW1hcnkocGNhKQpgYGAKCgpgYGB7cn0KYmlwbG90KHBjYSkKYGBgCgpgYGB7cn0KRnAgPC0gcGNhJHNjb3JlcwpHcyA8LSBwY2EkbG9hZGluZ3MKCkZzIDwtIEZwICUqJSBkaWFnKDEvcGNhJHNkZXYpCkdwIDwtIEdzICUqJSBkaWFnKHBjYSRzZGV2KSAqIDIKCmNvbC5jbGFzcyA8LSBhcy5udW1lcmljKHRyYWluJFY1OCkKY29sLmNsYXNzW2NvbC5jbGFzcz09MV0gPC0gInJlZCIKY29sLmNsYXNzW2NvbC5jbGFzcz09Ml0gPC0gImdyZWVuIgpjb2wuY2xhc3NbY29sLmNsYXNzPT0zXSA8LSAiYmx1ZSIKY29sLmNsYXNzW2NvbC5jbGFzcz09NF0gPC0gInllbGxvdyIKY29sLmNsYXNzW2NvbC5jbGFzcz09NV0gPC0gInB1cnBsZSIKCnBsb3QoRnNbLDFdLCBGc1ssMl0sIGFzcD0xLCBjb2wgPSBjb2wuY2xhc3MsIHhsYWIgPSAiRmlyc3QgcHJpbmNpcGFsIGNvbXBvbmVudCIsIHlsYWIgPSAiU2Vjb25kIHByaW5jaXBhbCBjb21wb25lbnQiKQphcnJvd3MocmVwKDAsZGltKEdzKVsxXSkscmVwKDAsZGltKEdzKVsxXSksIEdwWywxXSwgR3BbLDJdKQp0ZXh0KEdwWywxXSwgR3BbLDJdLCBuYW1lcyh0cmFpbl9udW0pLCBjb2wgPSAiYmxhY2siKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJncmVlbiIsICJibHVlIiwgInllbGxvdyIsICJwdXJwbGUiKSwgbGVnZW5kPWMoJzAnLCcxJywnMicsJzMnLCAnNCcpKQpgYGAKCgpWMSwgVjU5LCBWNTcgY3JlYXRlcyBtYW55IHByb2JsZW1zCgpgYGB7cn0KcHJvYmxlbWF0aWMgPC0gYygiVjU3IiwgIlY1OSIpCnRyYWluIDwtIHRyYWluWywgIShuYW1lcyh0cmFpbikgJWluJSBwcm9ibGVtYXRpYyldCmZkYSA8LSBsZGEoVjU4flYzK1Y0K1Y5K1YxMCtWMTErVjEyK1YxNCtWMTUrVjE2K1YxOCtWMTkrVjIzK1YyNCtWMjUrVjI2K1YyNytWMjkrVjMxK1YzMitWMzMrVjM0K1YzNStWMzcrVjM4K1YzOStWNDArVjQxK1Y0MytWNDQrVjUxK1Y1NitWNjArVjYxK1Y2MytWNjUrVjY3K1Y2OCwgZGF0YT10cmFpbikKI3Bsb3QoZmRhKQpsb2FkaW5ncyA8LSBwcmVkaWN0KGZkYSkkeApwbG90KGxvYWRpbmdzLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCmBgYAoKCmBgYHtyfQp0cmFpbiRMRDEgPC0gbG9hZGluZ3NbLDFdCnRyYWluJExEMiA8LSBsb2FkaW5nc1ssMl0KdHJhaW4kTEQzIDwtIGxvYWRpbmdzWywzXQp0cmFpbiRMRDQgPC0gbG9hZGluZ3NbLDRdCgpmZGFfdGVzdCA8LSBwcmVkaWN0KGZkYSwgbmV3ZGF0YSA9IHRlc3QpCnRlc3QkTEQxIDwtIGZkYV90ZXN0JHhbLDFdCnRlc3QkTEQyIDwtIGZkYV90ZXN0JHhbLDJdCnRlc3QkTEQzIDwtIGZkYV90ZXN0JHhbLDNdCnRlc3QkTEQ0IDwtIGZkYV90ZXN0JHhbLDRdCmBgYAoKCkFuw6FsaXNpcyBwb3IgY29ycmVzcG9uZGVuY2lhcwoKYGBge3J9CmFjIDwtIG1qY2EoY2xldlssbmFtZXMoY2xldikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iQnVydCIpCnBsb3QoYWMsIG1haW49Ik1DQSBiaXBsb3Qgb2YgQnVydCBtYXRyaXggd2l0aCBkYXRhIikKCmFjX2luZCA8LSBtamNhKHRyYWluWyxuYW1lcyh0cmFpbikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iaW5kaWNhdG9yIiwgcmV0aSA9IFQpCnBsb3QoYWNfaW5kJHJvd2Nvb3JkLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCgptY2EuZmVhdHVyZXMgPC0gYWNfaW5kJHJvd2Nvb3JkCmBgYAoKCiMjIyAzLiBSZXNhbXBsaW5nIHByb3RvY29sClByaW5jaXBhbCB1dGlsaXR5IGZ1bmN0aW9uIGZvciBkb2luZyBjcm9zcy12YWxpZGF0aW9uLiAKCkNyZWF0ZSB0aGUgQ1YgZnVuY3Rpb24gZm9yIGFueSBtb2RlbCB3aXRoIGFzc29jaWF0ZWQgcHJlZGljdCBhbmQgdXBkYXRlIGZ1bmN0aW9ucy4KCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uIDwtIGZ1bmN0aW9uKGRhdGEsIHRhcmdldCwgbW9kZWwsIHRpbWVzLCBuZm9sZHMsIG5lZWQuY2xhc3MpewogIHNldC5zZWVkKDAyMDIpCiAgQ1YuZm9sZHMgPC0gZ2VuZXJhdGVDVlJ1bnModGFyZ2V0LCBudGltZXM9dGltZXMsIG5mb2xkPW5mb2xkcywgc3RyYXRpZmllZD1UUlVFKQoKICBlcnIudG90YWwgPC0gYygpCiAgZm9yIChpIGluIDE6dGltZXMpewogICAgZXJyLm9uZXRpbWUgPC0gYygpCiAgICBmb3IgKGogaW4gMTpuZm9sZHMpewogICAgICBwcmludChwYXN0ZTAoIkZvbGQ6ICIsIGopKQogICAgICB2YWwgPC0gdW5saXN0KENWLmZvbGRzW1tpXV1bW2pdXSkKICAgICAgCiAgICAgIHRyIDwtIGRhdGFbLXZhbCxdCiAgICAgIHZhIDwtIGRhdGFbdmFsLF0KCiAgICAgIG1vZGVsIDwtIHVwZGF0ZShtb2RlbCwgZGF0YT10cikKICAgICAgcHJlZCA8LSBwcmVkaWN0KG1vZGVsLCBuZXdkYXRhPXZhKQogICAgICBpZiAobmVlZC5jbGFzcyl7CiAgICAgICAgcHJlZCA8LSBwcmVkJGNsYXNzCiAgICAgIH0KICAgICAgCiAgICAgIGVyci50YWJsZSA8LSB0YWJsZShUcnVlPXRhcmdldFt2YWxdLCBQcmVkaWN0ZWQ9cHJlZCkKICAgICAgZXJyIDwtIDEtc3VtKGRpYWcoZXJyLnRhYmxlKSkvc3VtKGVyci50YWJsZSkKICAgICAgZXJyLm9uZXRpbWUgPC0gYyhlcnIub25ldGltZSwgZXJyKQogICAgfQogICAgZXJyLnRvdGFsIDwtIGMoZXJyLnRvdGFsLCBtZWFuKGVyci5vbmV0aW1lKSkKICAgIHByaW50KHBhc3RlMCgiSXRlcmF0aW9uICIsIGksICIsIG1lYW4gZXJyb3I6ICIsIG1lYW4oZXJyLm9uZXRpbWUpKSkKICB9CiAgbWVhbihlcnIudG90YWwpCn0KYGBgCgoKQ3Jvc3MtdmFsaWRhdGlvbiBmb3Igay1OZWFyZXN0IE5laWdoYm91cgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24ua25uIDwtIGZ1bmN0aW9uKGRhdGEsIHRhcmdldCwgdGltZXMsIG5mb2xkcywgSyl7CiAgc2V0LnNlZWQoMDIwMikKICBDVi5mb2xkcyA8LSBnZW5lcmF0ZUNWUnVucyh0YXJnZXQsIG50aW1lcz10aW1lcywgbmZvbGQ9bmZvbGRzLCBzdHJhdGlmaWVkPVRSVUUpCgogIGVyci50b3RhbCA8LSBjKCkKICBmb3IgKGkgaW4gMTp0aW1lcyl7CiAgICBlcnIub25ldGltZSA8LSBjKCkKICAgIGZvciAoaiBpbiAxOm5mb2xkcyl7CiAgICAgIHZhbCA8LSB1bmxpc3QoQ1YuZm9sZHNbW2ldXVtbal1dKQogICAgICAKICAgICAgdHIgPC0gZGF0YVstdmFsLF0KICAgICAgdmEgPC0gZGF0YVt2YWwsXQoKICAgICAgcHJlZCA8LSBrbm4odHIsIHZhLCB0YXJnZXRbLXZhbF0sIGsgPSBLKQoKICAgICAgZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGFyZ2V0W3ZhbF0sIFByZWRpY3RlZD1wcmVkKQogICAgICBlcnIgPC0gMS1zdW0oZGlhZyhlcnIudGFibGUpKS9zdW0oZXJyLnRhYmxlKQogICAgICBlcnIub25ldGltZSA8LSBjKGVyci5vbmV0aW1lLCBlcnIpCiAgICB9CiAgICBlcnIudG90YWwgPC0gYyhlcnIudG90YWwsIG1lYW4oZXJyLm9uZXRpbWUpKQogIH0KICBtZWFuKGVyci50b3RhbCkKfQpgYGAKCgpDcm9zcy12YWxpZGF0aW9uIE5haXZlLUJheWVzCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbi5uYWl2ZSA8LSBmdW5jdGlvbihkYXRhLCB0YXJnZXQsIG1vZGVsLCB0aW1lcywgbmZvbGRzKXsKICBzZXQuc2VlZCgwMjAyKQogIENWLmZvbGRzIDwtIGdlbmVyYXRlQ1ZSdW5zKHRhcmdldCwgbnRpbWVzPXRpbWVzLCBuZm9sZD1uZm9sZHMsIHN0cmF0aWZpZWQ9VFJVRSkKCiAgZXJyLnRvdGFsIDwtIGMoKQogIGZvciAoaSBpbiAxOnRpbWVzKXsKICAgIGVyci5vbmV0aW1lIDwtIGMoKQogICAgZm9yIChqIGluIDE6bmZvbGRzKXsKICAgICAgcHJpbnQocGFzdGUwKCJGb2xkOiAiLCBqKSkKICAgICAgdmFsIDwtIHVubGlzdChDVi5mb2xkc1tbaV1dW1tqXV0pCiAgICAgIAogICAgICB0ciA8LSBkYXRhWy12YWwsXQogICAgICB2YSA8LSBkYXRhW3ZhbCxdCgogICAgICBtb2RlbCA8LSBuYWl2ZUJheWVzKFY1OH4uLUxEMS1MRDItTEQzLUxENCwgZGF0YT10cikKICAgICAgcHJlZCA8LSBwcmVkaWN0KG1vZGVsLCBuZXdkYXRhPXZhKQogICAgICAKICAgICAgZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGFyZ2V0W3ZhbF0sIFByZWRpY3RlZD1wcmVkKQogICAgICBlcnIgPC0gMS1zdW0oZGlhZyhlcnIudGFibGUpKS9zdW0oZXJyLnRhYmxlKQogICAgICBlcnIub25ldGltZSA8LSBjKGVyci5vbmV0aW1lLCBlcnIpCiAgICB9CiAgICBlcnIudG90YWwgPC0gYyhlcnIudG90YWwsIG1lYW4oZXJyLm9uZXRpbWUpKQogICAgcHJpbnQocGFzdGUwKCJJdGVyYXRpb24gIiwgaSwgIiwgbWVhbiBlcnJvcjogIiwgbWVhbihlcnIub25ldGltZSkpKQogIH0KICBtZWFuKGVyci50b3RhbCkKfQpgYGAKCgojIyMgNC4gTW9kZWxzClRoZSBtb2RlbHMgd2UgYXJlIGdvaW5nIHRvIHVzZSBhcmU6IAogIC0gTERBCiAgLSBRREEKICAtIFJEQQogIC0gay1OTgogIC0gTmHDr3ZlIEJheWVzCiAgLSBHTE0KIAoKYGBge3J9CnJkYS5tb2RlbCA8LSByZGEoVjU4flYzK1Y0K1Y5K1YxMCtWMTErVjEyK1YxNCtWMTUrVjE2K1YxOCtWMTkrVjIwK1YyMStWMjIrVjIzK1YyNCtWMjUrVjI2K1YyNytWMjkrVjMxK1YzMitWMzMrVjM0K1YzNStWMzcrVjM4K1YzOStWNDArVjQxK1Y0MytWNDQrVjUxK1Y1NitWNjArVjYxK1Y2MytWNjUrVjY3K1Y2OCwgZGF0YT10cmFpbikKbmFpdmUubW9kZWwgPC0gbmFpdmVCYXllcyhWNTh+VjMrVjQrVjkrVjEwK1YxMStWMTIrVjE0K1YxNStWMTYrVjE4K1YxOStWMjArVjIxK1YyMitWMjMrVjI0K1YyNStWMjYrVjI3K1YyOStWMzErVjMyK1YzMytWMzQrVjM1K1YzNytWMzgrVjM5K1Y0MCtWNDErVjQzK1Y0NCtWNTErVjU2K1Y2MCtWNjErVjYzK1Y2NStWNjcrVjY4LCBkYXRhPXRyYWluKQpgYGAKCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCByZGEubW9kZWwsIDEwLCAxMCwgVCkKYGBgCgoKYGBge3J9CnJkYS5tb2RlbC5mZGEgPC0gcmRhKFY1OH4uLGRhdGE9dHJhaW4pCmBgYAoKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIHJkYS5tb2RlbC5mZGEsIDEwLCAxMCwgVCkKYGBgCgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24ubmFpdmUodHJhaW4sIHRyYWluJFY1OCwgbmFpdmUubW9kZWwsIDEwLCAxMCkKYGBgCgoKYGBge3J9CmVyciA8LSBjKCkKZm9yIChrIGluIDE6MjApewogIGVyciA8LSBjKGVyciwgY3Jvc3MudmFsaWRhdGlvbi5rbm4odHJhaW4sIHRyYWluJFY1OCwgMTAsMTAsIGspKQp9CmBgYAoKCmBgYHtyfQpwbG90KGVyciwgdHlwZSA9ICJsIikKZXJyCmBgYAoKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uLmtubih0cmFpbiwgdHJhaW4kVjU4LCAxMCwgMTAsIDEpCgptdWx0aW5vbWlhbC5tb2RlbCA8LSBtdWx0aW5vbShWNTh+LiwgZGF0YT10cmFpbikKCmNyb3NzLnZhbGlkYXRpb24odHJhaW4sIHRyYWluJFY1OCwgbXVsdGlub21pYWwubW9kZWwsIDEwLCAxMCwgRikKCm11bHRpbm9taWFsLm1vZGVsLnN0ZXAgPC0gc3RlcChtdWx0aW5vbWlhbC5tb2RlbCkKCmNyb3NzLnZhbGlkYXRpb24odHJhaW4sIHRyYWluJFY1OCwgbXVsdGlub21pYWwubW9kZWwuc3RlcCwgMTAsIDEwLCBGKQoKbXVsdGlub21pYWwubW9kZWwubm9GREEgPC0gbXVsdGlub20oVjU4fi4tTEQxLUxEMi1MRDMtTEQ0LCBkYXRhPXRyYWluKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQSwgMTAsIDEwLCBGKQoKbXVsdGlub21pYWwubW9kZWwubm9GREEuc3RlcCA8LSBzdGVwKG11bHRpbm9taWFsLm1vZGVsLm5vRkRBKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQS5zdGVwLCAxMCwgMTAsIEYpCmBgYAoKCiMjIFRlc3QgZXJyb3IKCmBgYHtyfQpyZGEubW9kZWwgPC0gdXBkYXRlKHJkYS5tb2RlbC5mZGEsIGRhdGE9dHJhaW4pCnByZWQudGVzdCA8LSBwcmVkaWN0KHJkYS5tb2RlbC5mZGEsIHRlc3QpCnByZWQudGVzdCA8LSBwcmVkLnRlc3QkY2xhc3MKKGVyci50YWJsZSA8LSB0YWJsZShUcnVlPXRlc3QkVjU4LCBQcmVkPXByZWQudGVzdCkpCihlcnIudGVzdCA8LSAxLXN1bShkaWFnKGVyci50YWJsZSkpL3N1bShlcnIudGFibGUpKQpgYGAKCiMgUGFydGUgMgoKV2Ugd2lsbCBub3cgdXNlIG1vZGVsIG9uIGEgZGF0YXNldCBub3QgeWV0IGludmVzdGlnYXRlZDogSHVuZ2FyeS4KCmBgYHtyfQpodW5nIDwtIHJlYWQudGFibGUoIi4uL2RhdGEvcHJvY2Vzc2VkLmh1bmdhcmlhbi5kYXRhIiwgaGVhZGVyPUYsIHNlcD0nLCcsIG5hLnN0cmluZ3M9Ij8iKQpzdW1tYXJ5KGh1bmcpCmBgYAoKIyMgTWlzc2luZyB2YWx1ZXMgdHJlYXRtZW50CgpgYGB7cn0Kcm0udmFyIDwtIGMoIlYxMSIsICJWMTIiLCAiVjEzIikKaHVuZyA8LSBodW5nWywhKG5hbWVzKGh1bmcpICVpbiUgcm0udmFyKV0KYGBgCgpgYGB7cn0KaHVuZ1tpcy5uYShodW5nKV0gPC0gLTkKdmFyLmltcHV0YWJsZSA8LSBjKCJWNCIsICJWNSIsICJWNiIsICJWNyIsICJWOCIsICJWOSIpCmZvciAobm9tIGluIHZhci5pbXB1dGFibGUpewogIGh1bmdbLCBub21dW2h1bmdbLCBub21dID09IC05XSA8LSBOQQogIHZhcmlhYmxlIDwtIGh1bmdbLCBub21dCiAgaHVuZ1ssIG5vbV0gPC0ga25uLmltcHV0YXRpb24oaHVuZywgdmFyaWFibGUsIG5vbSwgNykKfQpgYGAKCmBgYHtyfQp2YXIuZmFjdG9yIDwtIGMoIlYyIiwgIlYzIiwgIlY2IiwgIlY3IiwgIlY5IiwgIlYxNCIpCmZvciAobm9tIGluIHZhci5mYWN0b3IpewogIGh1bmdbLCBub21dIDwtIGFzLmZhY3Rvcihhcy5jaGFyYWN0ZXIoaHVuZ1ssbm9tXSkpCn0KCiMgTGV2ZWwgMiBoYXMgZmV3IHZhbHVlcwpodW5nJFY3W2h1bmckVjcgPT0gMl0gPC0gMQpodW5nJFY3IDwtIGRyb3BsZXZlbHMoaHVuZyRWNykKbGV2ZWxzKGh1bmckVjcpCgojIExldmVsIDIgaGFzIG9uZSBhbm9tYWwgdmFsdWUKaHVuZyRWOVtodW5nJFY5ID09IDJdIDwtIDEKaHVuZyRWOSA8LSBkcm9wbGV2ZWxzKGh1bmckVjkpCmxldmVscyhodW5nJFY5KQpgYGAKCgpgYGB7cn0Kc3VtbWFyeShodW5nKQpgYGAKCiMjIEdhdXNzaWFuaXR5IHRyZWF0bWVudAoKYGBge3J9Cm5hbWVzX251bSA8LSBjKCkKZm9yKGkgaW4gMTpuY29sKGh1bmcpKXsKICBpZiAoIWlzLmZhY3RvcihodW5nWyxpXSkpIHsKICAgIG5hbWVzX251bSA8LSBjKG5hbWVzX251bSwgaSkKICB9Cn0KaHVuZ19udW1lcmljIDwtIGh1bmdbLG5hbWVzX251bV0KCmh1bmdfY29yIDwtIGNvcihodW5nX251bWVyaWMpCndoaWNoKGh1bmdfY29yID4gMC41ICYgaHVuZ19jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKd2hpY2goLWh1bmdfY29yID4gMC41ICYgLWh1bmdfY29yIDwgMSwgYXJyLmluZCA9IFRSVUUpCmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGh1bmdfbnVtZXJpYykpewogICAgcXFub3JtKGh1bmdfbnVtZXJpY1ssaV0sIG1haW4gPSBjKCJRLVEgUGxvdDogIiwgbmFtZXMoaHVuZ19udW1lcmljKVtpXSkpCiAgICBxcWxpbmUoaHVuZ19udW1lcmljWyxpXSwgY29sPTIpCn0KYGBgCgoKQm94Y294CgpgYGB7cn0KcGFyKG1mcm93ID0gYygyLDMpKQpmb3IoaSBpbiAxOmxlbmd0aChodW5nX251bWVyaWMpKXsKICBib3hjb3gobG0oaHVuZ19udW1lcmljWyxpXS1taW4oaHVuZ19udW1lcmljWyxpXSkrMWUtNn4xKSxsYW1iZGEgPSBzZXEoLTEsIDEuNSwgYnk9MC4xKSwgeGxhYiA9IGMobmFtZXMoaHVuZ19udW1lcmljKSlbaV0pCn0KYGBgCgpXZSBzZWUgdmFyIDEwIHJlcXVpcmVzIGEgbG9nYXJpdG1pYyB0cmFuc2Zvcm1hdGlvbi4KCmBgYHtyfQpodW5nWywiVjEwIl0gPC0gbG9nKDEraHVuZ1ssIlYxMCJdKQpgYGAgCgpgYGB7cn0KcXFub3JtKGh1bmdbLCJWMTAiXSkKcXFsaW5lKGh1bmdbLCJWMTAiXSwgY29sPTIpCmBgYAoKIyMgRmVhdHVyZSBleHRyYWN0aW9uCgpTZXBhcmUgdHJhaW4gYW5kIHRlc3QgZGF0YSwgc2VlZCBmb3IgcmVwcm9kdWNpYmlsaXR5LgoKYGBge3J9CnNldC5zZWVkKDIwMDApCm4gPC0gbnJvdyhodW5nKQp0cmFpbi5sZW5naHQgPC0gcm91bmQoMipuLzMpCgpodW5nIDwtIGh1bmdbc2FtcGxlKG4pLF0KdHJhaW4gPC0gaHVuZ1sxOnRyYWluLmxlbmdodCxdCnRlc3QgPC0gaHVuZ1sodHJhaW4ubGVuZ2h0KzEpOm4sXQpgYGAKCmBgYHtyfQpuYW1lc19udW0gPC0gYygpCmZvcihpIGluIDE6bmNvbCh0cmFpbikpewogIGlmICghaXMuZmFjdG9yKHRyYWluWyxpXSkpIHsKICAgIG5hbWVzX251bSA8LSBjKG5hbWVzX251bSwgaSkKICB9Cn0KdHJhaW5fbnVtIDwtIHRyYWluWyxuYW1lc19udW1dCmBgYAoKRXh0cmFjdCBQQ0EgZmVhdHVyZXMuCgpgYGB7cn0KcGNhIDwtIHByaW5jb21wKHRyYWluX251bSkKc2NyZWVwbG90KHBjYSkKc3VtbWFyeShwY2EpCmBgYAoKCmBgYHtyfQpiaXBsb3QocGNhKQpgYGAKCmBgYHtyfQpGcCA8LSBwY2Ekc2NvcmVzCkdzIDwtIHBjYSRsb2FkaW5ncwoKRnMgPC0gRnAgJSolIGRpYWcoMS9wY2Ekc2RldikKR3AgPC0gR3MgJSolIGRpYWcocGNhJHNkZXYpICowLjMKCmNvbC5jbGFzcyA8LSBhcy5udW1lcmljKHRyYWluJFYxNCkKY29sLmNsYXNzW2NvbC5jbGFzcz09MF0gPC0gInJlZCIKY29sLmNsYXNzW2NvbC5jbGFzcz09MV0gPC0gImJsdWUiCgpwbG90KEZzWywxXSwgRnNbLDJdLCBhc3A9MSwgY29sID0gY29sLmNsYXNzLCB4bGFiID0gIkZpcnN0IHByaW5jaXBhbCBjb21wb25lbnQiLCB5bGFiID0gIlNlY29uZCBwcmluY2lwYWwgY29tcG9uZW50IikKYXJyb3dzKHJlcCgwLGRpbShHcylbMV0pLHJlcCgwLGRpbShHcylbMV0pLCBHcFssMV0sIEdwWywyXSkKdGV4dChHcFssMV0sIEdwWywyXSwgbmFtZXModHJhaW5fbnVtKSwgY29sID0gImJsYWNrIikKbGVnZW5kKCJib3R0b21yaWdodCIsIGZpbGw9YygicmVkIiwiYmx1ZSIpLCBsZWdlbmQ9YygnMCcsJzEnKSkKYGBgCgpgYGB7cn0KZmRhIDwtIGxkYShWMTR+LiwgZGF0YT10cmFpbikKI3Bsb3QoZmRhKQpsb2FkaW5ncyA8LSBwcmVkaWN0KGZkYSkkeApwbG90KGxvYWRpbmdzLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImJsdWUiKSwgbGVnZW5kPWMoJzAnLCcxJykpCmBgYAoKCmBgYHtyfQp0cmFpbiRMRDEgPC0gbG9hZGluZ3NbLDFdCgpmZGFfdGVzdCA8LSBwcmVkaWN0KGZkYSwgbmV3ZGF0YSA9IHRlc3QpCnRlc3QkTEQxIDwtIGZkYV90ZXN0JHhbLDFdCmBgYAoKCkFuw6FsaXNpcyBwb3IgY29ycmVzcG9uZGVuY2lhcwoKYGBge3J9CmFjIDwtIG1qY2EoaHVuZ1ssbmFtZXMoaHVuZykgJWluJSB2YXIuZmFjdG9yXSwgbGFtYmRhPSJCdXJ0IikKcGxvdChhYywgbWFpbj0iTUNBIGJpcGxvdCBvZiBCdXJ0IG1hdHJpeCB3aXRoIGRhdGEiKQoKYWNfaW5kIDwtIG1qY2EodHJhaW5bLG5hbWVzKHRyYWluKSAlaW4lIHZhci5mYWN0b3JdLCBsYW1iZGE9ImluZGljYXRvciIsIHJldGkgPSBUKQpwbG90KGFjX2luZCRyb3djb29yZCwgY29sID0gY29sLmNsYXNzKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJibHVlIiksIGxlZ2VuZD1jKCcwJywnMScpKQoKbWNhLmZlYXR1cmVzIDwtIGFjX2luZCRyb3djb29yZApgYGAKCiMjIAo=